// Copyright 2022 Jeroen van Nugteren

// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
// OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#ifndef DM_DISTMESH2D_HH
#define DM_DISTMESH2D_HH

// general headers
#include <cassert>
#include <memory>
#include <json/json.h>
#include <armadillo> 

// fx-common headers
#include "rat/common/extra.hh"
#include "rat/common/error.hh"

// distmesh-cpp headers
#include "rat/common/gmshfile.hh"
#include "rat/common/node.hh"
#include "dfones.hh"
#include "distfun.hh"


// code specific to Rat
namespace rat{namespace dm{

	// shared pointer definition
	typedef std::shared_ptr<class DistMesh2D> ShDistMesh2DPr;

	// a lightweight tri and quad mesh generator
	// based on the 2D distmesh from Per Olof Persson
	// http://persson.berkeley.edu/distmesh/
	// Per-Olof Persson and Gilbert Strang, "A Simple Mesh 
	// Generator in MATLAB," SIAM Review Vol. 46 (2) 2004.
	// using delaunay triangulation algorithm from
	// Paul Bourke http://paulbourke.net/papers/triangulate/
	// both algorithms have been translated to Armadillo style
	// programming
	class DistMesh2D: public cmn::Node{
		// input parameters
		protected:
			// tri-mesh or quad mesh
			bool quad_ = true; // default when used in Rat

			// element size settings
			fltp h0_;

			// fixed node positions [x,y]
			arma::Mat<fltp> pfix_ext_ = arma::Mat<fltp>(0,2);

			// distance function
			// negative inside positive outside
			ShDistFunPr fd_;

			// scaling function
			ShDistFunPr fh_;

			// soft limit on the number of nodes
			// this can safely be changed it is just 
			// to prevent you from going out of memory
			// when forgetting to set the element size
			arma::uword num_node_limit_ = 500000;

			// seed for rng
			unsigned int rng_seed_ = 1001;
			fltp rng_strength_ = RAT_CONST(0.2); // applies random offset to the starting positions of the vertices

		// mesher settings
		// do not touch
		protected:
			// tolerances
			fltp dptol_ = RAT_CONST(0.001); 
			fltp ttol_ = RAT_CONST(0.1); 
			fltp Fscale_ = RAT_CONST(1.2); 
			fltp deltat_ = RAT_CONST(0.2); 

			// quad recombination
			fltp recombination_treshold_ = RAT_CONST(0.1);
			bool alternate_recombination_ = false;

			// iterations
			arma::uword densityctrlfreq_ = 30;
			arma::uword maxiter_ = 1000;

		// internal storage
		protected:
			// derived tolerances
			fltp deps_; // = std::sqrt(arma::datum::eps)*h0
			fltp geps_; // = 0.001*h0; 

			// fixed points
			arma::Mat<fltp> pfix_;

			// bounding box [x,y]
			arma::Mat<fltp>::fixed<2,2> bbox_;

			// stored mesh
			arma::Mat<arma::uword> qt_; // quads or triangles [p1,p2,p3,(p4)]
			arma::Mat<fltp> p_; // 2D points [x,y]

			// edges and surface edges
			arma::Mat<arma::uword> e_;
			arma::Mat<arma::uword> es_; 
			arma::Mat<arma::uword> ei_; 

			// surface nodes
			arma::Col<arma::uword> ps_;


		// methods
		public:
			// constructcor
			DistMesh2D();

			// destructor
			virtual ~DistMesh2D(){};

			// factory
			static ShDistMesh2DPr create();


			// getting and setting
			// set element size
			void set_h0(const fltp h0);

			// set seed and strength for random number generator
			void set_rng_seed(const unsigned int rng_seed);
			void set_rng_strength(const fltp rng_strength);

			// recombine triangles to quad
			void set_quad(const bool quad);

			// set distance function
			void set_distfun(const ShDistFunPr &fd);

			// set scaling function
			void set_scalefun(const ShDistFunPr &fh);

			// set maximum number of iterations
			void set_maxiter(const arma::uword maxiter);

			// set fixed points
			void set_fixed(const arma::Mat<fltp> &pfix_ext);

			// set number of elements limit
			void set_num_node_limit(const arma::uword num_node_limit);

			// get mesh, nodes and elements
			const arma::Mat<arma::uword>& get_elements() const;
			const arma::Mat<fltp>& get_nodes() const;
			const arma::Mat<arma::uword>& get_surface_edges() const;
			const arma::Mat<arma::uword>& get_internal_edges() const;
			const arma::Col<arma::uword>& get_surface_node_idx() const;

			// get functions
			const ShDistFunPr& get_distfun()const;
			const ShDistFunPr& get_scalefun()const;


			// main mesher
			// create edge list
			void create_edges();

			// algorithm for converting a tri-mesh to a quad only mesh
			void tri2quad();

			// algorithm for removing high valence nodes in quad mesh
			void high_valence_removal(const arma::uword valence);

			// basic extrusion algorithm for making 3D meshes
			void extrude(const fltp ell, const arma::uword num_z);

			// smoothing (used for quad mesh only)
			void smoothen();

			// method for removing duplicates from a mesh
			static void remove_unreferenced(arma::Mat<fltp>&p, arma::Mat<arma::uword>&qt);
			void remove_duplicate();
			void fix_clockwise();

			// main distmesh algorithm
			void setup();

			// export to gmsh
			void export_gmsh(const cmn::ShGmshFilePr& gmsh) const;
			void export_gmsh_perimeter(const cmn::ShGmshFilePr& gmsh)const;

			// functionf for calculating area
			fltp calc_area() const;

			// get valence for each node
			arma::Col<arma::uword> get_node_valence()const;
			
			

			// delaunay functions
			// function for checking points in circle
			static bool circumcircle(
				arma::Row<fltp>::fixed<2> &pc, 
				fltp &r, 
				const arma::Row<fltp>::fixed<2> &p, 
				const arma::Mat<fltp>::fixed<3,2> &ptri);

			// delaunay triangulation function
			static arma::Mat<arma::uword> triangulation(
				const arma::Mat<fltp> &p);

			// section find math function
			static arma::Mat<arma::uword> find_sections(
				const arma::Row<arma::uword> &M);


			// distmesh helper functions
			// meshgrid for generating initial coordinates
			static void meshgrid(
				arma::Mat<fltp> &x, 
				arma::Mat<fltp> &y, 
				const arma::Row<fltp> &xgv, 
				const arma::Col<fltp> &ygv);

			// triangular to quad helper functions		
			// function for creating graph of all neighbouring elements
			static void element_graph(
				arma::Col<arma::uword> &row, 
				arma::Col<arma::uword> &col, 
				const arma::Mat<arma::uword> &t);

			// function for calculating the quality of elements
			static arma::Col<fltp> element2quality(
				const arma::Mat<arma::uword> &tq, 
				const arma::Mat<fltp> &p);
			static arma::Col<fltp> element2quality_alt(
				const arma::Mat<arma::uword> &tq, 
				const arma::Mat<fltp> &p);

			// algorithm for merging two triangles into one quad
			static arma::Row<arma::uword>::fixed<4> merge2triangles(
				const arma::Row<arma::uword>::fixed<3> &t1, 
				const arma::Row<arma::uword>::fixed<3> &t2);

			// algorithm for performing catmull clark subdivision
			static void refine(
				arma::Mat<arma::uword> &qr, 
				arma::Mat<fltp> &pr, 
				const arma::Mat<arma::uword> &qt, 
				const arma::Mat<fltp> &p);


			// serialization and deserializetion of input settings
			static std::string get_type();
			void serialize(
				Json::Value &js, 
				cmn::SList &list) const override;
			void deserialize(
				const Json::Value &js, 
				cmn::DSList &list, 
				const cmn::NodeFactoryMap &factory_list, 
				const boost::filesystem::path &pth) override;
	};

}}

#endif
